method_mean = 'riemann';
method_dist = 'riemann';

%% segement all the files
FileList = {
   
  %  'mohEEG111.mat'
      'Marion12_2504.mat'
    %  'Tim12_2804.mat'
    %  'Romain12_2904.mat'
    };

%FileList = { 'session2.acq.mat' };
%Fs = 2000; default
Fs = 2000;

window = 6000;
window_movement = 6000;
cov_mat_per_class = 10;%number of covariance matrices per movement
% also this is the number of overlaping windows per movement 
windows_sizes = floor( window_movement / cov_mat_per_class );
overlap = window_movement - window_movement / 8;
shift = 0;
offset_train = 500;

X = []; % contains the trials from all classes
Y = []; % contains the classification from all classes

Nclass = 12;
CallClass = cell(Nclass,cov_mat_per_class);
n_channels = 10;

%start_of_class = [];
%end_of_class = [];


for curFileIx = 1:length(FileList)
   


        load(['D:\temp\subjects\' FileList{curFileIx}]);
        disp(['Subject file:' FileList{curFileIx}]);

        datar = data(:,1:size(data,2)-2); % remove 2 trigger channel at the end

        %trial_channel = data(:,8);
        trial_channel = data(:, size(data,2)-1);

        %class_separation_channel = data(:,9);
        class_separation_channel = data(:, size(data,2));

        classN = length(count_sequence(class_separation_channel,5));  % "5" tags each class


        classRest = classN; % the rest class is always the last

        trialsN = length(count_sequence(trial_channel,5));  % "5" tags each trial
            
        disp(['Classes/Movements detected: ' int2str(classN)]);
        disp(['Total trails detected: ' int2str(trialsN)]);
        disp(['Frequency: ' int2str(Fs)]);

        triggers_start_per_class = zeros(Nclass,20);
        
        for m = 1:classN
            
            %Pulsar BIOPAC 
            %[best_index,best_trigger] = generate_trials_postitions_using_discrete_channel(m, trial_channel,class_separation_channel);
            [best_index,best_trigger1,class_start_pos,class_end_pos] = generate_trials_postitions_using_discrete_channel(m, trial_channel,class_separation_channel);
            %start_of_class = [start_of_class class_start_pos];
            %end_of_class = [end_of_class class_end_pos];
               
            % cut signal
            %index = find(diff(trigger)==1); % because we do it earlier now!!!!
            trials_per_class = zeros(size(datar,2),window,length(best_index));
            disp(['Trials: ' int2str(length(best_index)) ' ' FileList{curFileIx}])

            %get the data for these triggers
            for i=1:length(best_index)
                %calculate a range around the trigger
                buffer = datar(best_index(i)+shift:best_index(i)+shift+window-1,:);
                trials_per_class(:,:,i) = buffer';
                
                triggers_start_per_class(m,i) = best_index(i)+shift;
            end

            X = cat(3,X,trials_per_class); % put all trials together
            Y = cat(1,Y,ones(length(best_index),1) * m);

        end;
        
        %all data is available here for the subject for each class
        
        %COVtest = zeros(size(datar,2),size(datar,2),size(X,3),cov_mat_per_class);
    disp('--------------------------------');% next file
end

%% Perform training - calcualte barycenters
for c = 1:Nclass

    epochs_per_class = X(:,:,Y==c);

    covv = zeros(size(datar,2),size(datar,2),size(epochs_per_class,3));            

    for j=1:cov_mat_per_class % cut each epoch on several consequal chunks of signal

         for i = 1:size(epochs_per_class,3) % go through all signal epochs


            %calculate cov for the this part of the signal
            buffer = epochs_per_class(:, (j-1) * windows_sizes + 1 : j * windows_sizes ,i);
            covv(:,:,i) = covariances(buffer);

            %calculate the barycenter

         end;

         CallClass{c,j} = mean_covariances(covv,method_mean);
     end;

end;

%% Perform linear training - divide the window and calculate the distances for each part 

[epochs epochsIndex] = eeg2epochs_index(datar', window_movement ,overlap);

NTesttrial = size(epochs,3);

XX = zeros(NTesttrial, Nclass * cov_mat_per_class);
%yy = zeros(NTesttrial,1); 
yy = ones(NTesttrial,1) * 13; % default class = 13

taken_as_movement= 0;

for i=1:NTesttrial

            for c=1:Nclass
               
                 for j=1:cov_mat_per_class
                 
                    buffer = epochs(:, (j-1) * windows_sizes + 1 : j * windows_sizes ,i);
                    covv = covariances(buffer);
                    XX(i,(c-1) * cov_mat_per_class + j) = distance(covv, CallClass{c,j}, method_dist);
                    
                 end;
                 
                 for k=1:length(best_index)
                
                    if abs( (epochsIndex(i) + window_movement) - (triggers_start_per_class(c,k) + window_movement)) < offset_train
                       yy(i) = c;
                       taken_as_movement = taken_as_movement + 1;
                    end;
                 end
            end;         
end

featTr = XX;
Y2 = yy;
%featTe = XX;

%% Actual linear traininig based on the previously calculated distances for each window
lin = train(Y2,sparse(featTr),'-s 0 -q');

disp('Training completed');

%% Perform online

FileList = {
   
  %  'mohEEG111.mat'
      'Marion12_2504.mat'
     % 'Tim12_2804.mat'
     % 'Romain12_2904.mat'
    };


%load a new subject 
for curFileIx = 1:length(FileList)
   
        load(['D:\temp\subjects\' FileList{curFileIx}]);
        disp(['Subject file:' FileList{curFileIx}]);

        datar = data(:,1:size(data,2)-2); % remove 2 trigger channel at the end

        %trial_channel = data(:,8);
        trial_channel = data(:, size(data,2)-1);

        %class_separation_channel = data(:,9);
        class_separation_channel = data(:, size(data,2));

        classN = length(count_sequence(class_separation_channel,5));  % "5" tags each class

        classRest = classN; % the rest class is always the last

        trialsN = length(count_sequence(trial_channel,5));  % "5" tags each trial
            
        disp(['Classes/Movements detected: ' int2str(classN)]);
        disp(['Total trails detected: ' int2str(trialsN)]);
        disp(['Frequency: ' int2str(Fs)]);

        triggers_start_per_class = zeros(Nclass,20);
        
        for m = 1:classN
            
            %Pulsar BIOPAC 
            %[best_index,best_trigger] = generate_trials_postitions_using_discrete_channel(m, trial_channel,class_separation_channel);
            [best_index,best_trigger1,class_start_pos,class_end_pos] = generate_trials_postitions_using_discrete_channel(m, trial_channel,class_separation_channel);
            %start_of_class = [start_of_class class_start_pos];
            %end_of_class = [end_of_class class_end_pos];
               
            % cut signal
            %index = find(diff(trigger)==1); % because we do it earlier now!!!!
            trials_per_class = zeros(size(datar,2),window,length(best_index));
            disp(['Trials: ' int2str(length(best_index)) ' ' FileList{curFileIx}])

            %get the data for these triggers
            for i=1:length(best_index)            
                triggers_start_per_class(m,i) = best_index(i)+shift;
            end

        end;
    disp('--------------------------------');% next file
end
        
%produce epochs
[epochs epochsIndex] = eeg2epochs_index(datar', window_movement ,overlap);

NTesttrial = size(epochs,3);

XX = zeros(NTesttrial, Nclass * cov_mat_per_class);
%yy = zeros(NTesttrial,1); 
yy = ones(NTesttrial,1) * 13; % default class = 13

taken_as_movement= 0;

%calculate XX for new subject
for i=1:NTesttrial

            for c=1:Nclass
               
                 for j=1:cov_mat_per_class
                 
                    buffer = epochs(:, (j-1) * windows_sizes + 1 : j * windows_sizes ,i);
                    covv = covariances(buffer);
                    XX(i,(c-1) * cov_mat_per_class + j) = distance(covv, CallClass{c,j}, method_dist);
                    
                 end;
                 
                 for k=1:length(best_index)
                
                    if abs( (epochsIndex(i) + window_movement) - (triggers_start_per_class(c,k) + window_movement)) < offset_train
                       yy(i) = c;
                       taken_as_movement = taken_as_movement + 1;
                    end;
                 end
            end;         
end

%apply XX for the new subject


%featTr = XX;
Y2 = yy;
featTe = XX;

out1 = lin.w*featTe' + lin.bias;

[~,out1] = max(out1);
   
yte = lin.Label(out1);

disp('Online training completed');

correct = sum((yy == yte));
correct_mv = 0;

for i = 1:length(yy)
    if (yte(i)~=13 && yte(i)~=0 && yte(i) == yy (i) )
        correct_mv = correct_mv + 1;
    end;
end;

acc_movement = correct_mv / taken_as_movement
acc = correct / length(yy)

figure;
ytte = yte;
ytte(yte==13) = 0; % only for display
plot(ytte);
taken_as_movement
%[predicted_label] = predict(testing_label_vector, sparse(XX), lin, 'liblinear_options','col')
% with Y2 the labels
% FeatTr the training set (XX)
% and FeatTe the test set

